Chasing Peaks

# Libraries
library(tidyverse)
library(gmRi)
library(scales)
library(gt)
library(patchwork)
library(gtsummary)
library(gt)
library(sizeSpectra)
library(gganimate)
library(glue)
library(merTools)
conflicted::conflict_prefer("select", "dplyr")
conflicted::conflict_prefer("filter", "dplyr")

# Processing functions
source(here::here("Code/R/Processing_Functions.R"))

# ggplot theme
theme_set(
  theme_gmri(
    axis.line.y = element_line(),
    axis.ticks.y = element_line(), 
    rect = element_rect(fill = "white", color = NA),
    panel.grid.major.y = element_blank(),
    strip.text.y = element_text(angle  = 0),
    axis.text.x = element_text(size = 8),
    axis.text.y = element_text(size = 8),
    legend.position = "bottom"))



# labels for factor levels
area_levels <- c("GoM", "GB", "SNE", "MAB")
area_levels_long <- c("Gulf of Maine", "Georges Bank", "Southern New England", "Mid-Atlantic Bight")
area_levels_all <- c("Northeast Shelf", area_levels)
area_levels_long_all <- c("Northeast Shelf", area_levels_long)

# table to join for swapping shorthand for long-hand names
area_df <- data.frame(
  area = c("Scotian Shelf", "Gulf of Maine", "Georges Bank", "Southern New England", "Mid-Atlantic Bight", "Northeast Shelf"),
  survey_area = c("SS", "GoM", "GB", "SNE", "MAB", "Northeast Shelf"),
  area_titles = c("Scotian Shelf", "Gulf of Maine", "Georges Bank", "Southern New England", "Mid-Atlantic Bight", "Northeast Shelf"))


# Degree symbol
deg_c <- "\u00b0C"

Determining Spectra Minimum Sizes

Size spectra typically focus on the descending slope of the size distribution. This doc covers how that minima is selected and the process for doing so, and then the consequence of shifting the minima of the bounded distribution.

# Data used for Wigley estimates
wigley_in <- read_csv(here::here("Data/processed/wigley_species_trawl_data.csv")) %>% left_join(area_df)

The general idea is to bin the individual abundances by log2 size-bins, then determine the peak.

The following functions will aid in the labeling of log2 bins:

# Broad Distribution
#log2 bins for easy-of-access
#' @title Build Log 2 Bin Structure Dataframe
#' 
#' @description Used to build a dataframe containing equally spaced log2 bins for
#' size spectra analysis. Contains details on the left and right limits, midpoint, bin width, 
#' and a text label for the bins. log2bin number ascends with increasing size for eeasy plotting.
#'
#' @param log2_min 
#' @param log2_max 
#' @param log2_increment 
#'
#' @return
#' @export
#'
#' @examples
define_log2_bins <- function(log2_min = 0, log2_max = 13, log2_increment = 1){
  
  # How many bins
  n_bins  <- length(seq(log2_max, log2_min + log2_increment, by = -log2_increment))
  
  # Build Equally spaced log2 bin df
  log2_bin_structure <- data.frame(
    "log2_bins" = as.character(seq(n_bins, 1, by = -1)),
    "left_lim"  = seq(log2_max - log2_increment, log2_min, by = -log2_increment),
    "right_lim" = seq(log2_max, log2_min + log2_increment, by = -log2_increment)) %>% 
    mutate(
      bin_label    = str_c(round(2^left_lim, 3), " - ", round(2^right_lim, 3), "g"),
      bin_width    = 2^right_lim - 2^left_lim,
      bin_midpoint = (2^right_lim + 2^left_lim) / 2) %>% 
    arrange(left_lim)
  
  return(log2_bin_structure)
}






#' @title Assign Manual log2 Bodymass Bins - By weight
#'
#' @description Manually assign log2 bins based on individual length-weight bodymass 
#' in increments of 1 on the log2 scale. Returns data with bins assigned based on individual
#' length-weight biomass
#' 
#' Uses maximum weight, and its corresponding bin as the limit.
#'
#' @param wmin_grams Catch data prepared for mle calculation, use prep_wmin_wmax
#' @param log2_increment Equally spaced increments to use for log 2 bin sizes. Default = 0.5.
#'
#' @return
#' @export
#'
#' @examples
assign_log2_bins <- function(wmin_grams, log2_increment = 1){
  
  
  #### 1. Set up bodymass bins
  
  # filter missing weights
  size_data <- wmin_grams %>% 
    filter(wmin_g > 0,
           is.na(wmin_g) == FALSE,
           wmax_g > 0,
           is.na(wmax_g) == FALSE)
  
  # Get bodymass on log2() scale
  size_data$log2_weight <- log2(size_data$ind_weight_g)
  
  # Set up the bins - Pull min and max weights from data available
  #min_bin <- floor(min(size_data$log2_weight))
  min_bin <- 0
  max_bin <- ceiling(max(size_data$log2_weight))
  
  
  # Build a bin key, could be used to clean up the incremental assignment or for apply style functions
  log2_bin_structure <- define_log2_bins(
    log2_min = min_bin, 
    log2_max = max_bin, 
    log2_increment = log2_increment)
  
  
  
  # Loop through bins to assign the bin details to the data
  log2_assigned <- log2_bin_structure %>%
    split(.$log2_bins) %>%
    map_dfr(function(log2_bin){
      
      # limits and labels
      l_lim   <- log2_bin$left_lim
      r_lim   <- log2_bin$right_lim
      bin_num <- as.character(log2_bin$log2_bin)
      
      # assign the label to the appropriate bodymasses
      size_data %>% mutate(
        log2_bins = ifelse( between(log2_weight, l_lim, r_lim), bin_num, NA),
        log2_bins = as.character(log2_bins)) %>%
        drop_na(log2_bins)
      
    })
  
  # Join in the size bins
  log2_assigned <- left_join(log2_assigned, log2_bin_structure, by = "log2_bins")
  
  # return the data with the bins assigned
  return(log2_assigned)
  
}








#' @title Calculate Normalized and De-Normalized Abundances
#'
#' @description For binned size spectra estimation we use the stratified abundance divided by the
#' bin widths (normalized size spectra). Another way to present the data is to de-normalize, which 
#' takes those values and multiplies them by the mid-point of the log-scale bins.
#' 
#' min/max & bin_increments are used to enforce the presence of a size bin in the event that 
#' there is no abundance. This is done for comparing across different groups/areas that should 
#' conceivably have the same size range sampled.
#'
#' @param log2_assigned size data containing the bin assignments to use
#' @param min_log2_bin Minimum 2^x value for the size spectra being measured (>=)
#' @param max_log2_bin Maximum 2^x value for the size spectra being measured (<)
#' @param bin_increment The bin-width on log scale that separates each bin
#' @param ... Additional grouping factors with which to aggregate on besides the size bins themselves
#'
#' @return
#' @export
#'
#' @examples
aggregate_log2_bins <- function(
    log2_assigned, 
    min_log2_bin = 0, 
    max_log2_bin = 13, 
    bin_increment = 1,
    ...){
  
  # Full Possible Bin Structure
  # Fills in any gaps
  log2_bin_structure <- define_log2_bins(
    log2_min       = min_log2_bin, 
    log2_max       = max_log2_bin, 
    log2_increment = bin_increment)
  
  
  # Capture all the group levels with a cheeky expand()
  if(missing(...) == FALSE){
    log2_bin_structure <- log2_bin_structure %>% 
      tidyr::expand(left_lim, distinct(log2_assigned, ...)) %>% 
      full_join(log2_bin_structure)
  }
  
  
  
  # Get bin breaks
  log2_breaks <- sort(
    unique(c(log2_bin_structure$left_lim, log2_bin_structure$right_lim)))
  
  
  # Get Totals for bodymass and abundances
  log2_aggregates <- log2_assigned %>% 
    group_by(left_lim, ...) %>% 
    summarise(abundance   = sum(numlen_adj, na.rm = T),
              weight_g    = sum(wmin_g, na.rm = T),
              .groups = "drop")
  
  
  # Join back in what the limits and labels are
  # The defined bins and their labels enforce the size limits
  log2_prepped <- left_join(
    x = log2_bin_structure, 
    y = log2_aggregates)
  
  
  #### Fill Gaps with Zero's?? 
  # This ensures that any size bin that is intended to be in use is actually used
  log2_prepped <- log2_prepped %>% 
    mutate(
      abundance = ifelse(is.na(abundance), 1, abundance),
      weight_g = ifelse(is.na(weight_g), 1, weight_g))
  
  
  #### normalize abundances using the bin widths
  log2_prepped <- log2_prepped %>% 
    mutate(
      normalized_abund = abundance / bin_width,
      normalized_biom = weight_g / bin_width,
      # # Remove Bins Where Normalized Biomass < 0? No!
      # normalized_abund = ifelse(normalized_abund < 2^0, NA, normalized_abund),
      # norm_strat_abund = ifelse(norm_strat_abund < 2^0, NA, norm_strat_abund)
    )
  
  # Add de-normalized abundances (abundance * bin midpoint)
  log2_prepped <- log2_prepped %>% 
    mutate(
      denorm_abund = normalized_abund * bin_midpoint,
      denorm_biom = normalized_biom * bin_midpoint)
  
  # Return the aggregations
  return(log2_prepped)
  
}

Using those functions we can assign the bins based on the body-size column in the dataframe we would feed to our spectra slope workflow.

# Assign l2 bins
wig_l2 <- assign_log2_bins(wigley_in, log2_increment = 1)

# Create Shelf aggregate
wig_l2 <- wig_l2 %>% 
  mutate(area = "Northeast Shelf", 
         survey_area = "Northeast Shelf") %>% 
  bind_rows(wig_l2) %>% 
  mutate(area = factor(area, levels = area_levels_long_all),
         season = factor(season, levels = c("Spring", "Fall")))

From there, we can use one group as an example. In this plot peak normalized abundance is at the 2-4g bin.

This is not the most straightforward example because its multi-modal, but this will happen so it is instructive to use.

This example also highlights the common situation where normalized abundance approaches 0, which I’m still unsure how much that matters.

# Group details
glabs <- list(
  "yr" = 2000,
  "area" = "MAB",
  "season" = "Spring"
)

# Filter that out
test_case <- wig_l2 %>% 
  filter(
    est_year == glabs$yr, 
    survey_area == glabs$area, 
    season == glabs$season)



# Plot the thing - cheater way
ggplot(test_case) +
  geom_histogram(
    aes(
      x = wmin_g, 
      weight = (numlen_adj/bin_width),
      color = after_stat(count) == max(after_stat(count))), 
    breaks = 2^c(0:15), 
    linewidth = 1) +
  scale_color_manual(values = c("white", "orange")) +
  scale_y_continuous(
    trans = transform_log(base = 2),
    labels = label_log(base = 2),
    breaks = 2^seq(-12,10, 2)) +
  scale_x_continuous(
    trans = transform_log(base = 2),
    labels = label_log(base = 2),
    breaks = 2^c(0:15)) +
  labs(
    title = "Group Details:", 
    subtitle = str_c("Area: ", glabs$area,"\nSeason: ", glabs$season, "\nYear: ", glabs$yr ),
    y = "Normalized Abundance",
    x = "Bodyweight (g)")

I can’t auto-locate the tallest bin within all of that geom_histogram plot code, so it will be easier to do it outside and animate them using geom_col.

Here are all of them as an animation:

# Get the plotting summaries
wig_l2_aggs <- wig_l2 %>%
  unite("groupvar", area, season, est_year, sep = "XX") %>%
  split(.$groupvar) %>%
  map_dfr(function(x){

    # Get aggregates
    l2_aggs <- aggregate_log2_bins(x, bin_increment = 1)

    # Locate Tallest
    tallest <- l2_aggs %>%
      arrange(desc(normalized_abund)) %>%
      slice(1) %>%
      pull(left_lim)

    # Add info back to aggregates
    l2_aggs <- l2_aggs %>%
      mutate(is_tallest = if_else(left_lim == tallest, T, F))

    # Return
    return(l2_aggs)}, .id = "groupvar")  %>%
  separate(groupvar, into = c("area", "season", "est_year"), sep = "XX") %>%
  mutate(
    area = factor(area, levels = area_levels_long_all),
    season = factor(season, levels = c("Spring", "Fall")))






# # ggplot with animation
# animated_plot <- ggplot(wig_l2_aggs) +
#   geom_col(
#     aes(
#       x = left_lim,
#       y = normalized_abund,
#       fill = season,
#       color = is_tallest),
#     alpha = 0.6,
#     linewidth = 1.5) +
#   scale_y_continuous(
#     trans = transform_log(base = 2),
#     labels = label_log(base = 2),
#     breaks = 2^seq(-10,16, 2)) +
#   scale_x_continuous(
#     labels = label_math(expr = 2^.x),
#     breaks = c(0:15)) +
#   scale_fill_gmri() +
#   scale_color_manual(values = c("transparent", "black")) +
#   facet_grid(area ~ season, scale = "free") +
#   labs(
#     y = "Normalized Abundance",
#     x = "Bodyweight (g)",
#     title = 'Year: {closest_state}', # Dynamic title for year
#     fill = "Season") +
#   transition_states(
#     est_year,
#     transition_length = 1,
#     state_length = 5) +
#   ease_aes('linear')  # Smooth transition
# 
# 
# 
# # Animate and save
# animate(animated_plot, nframes = 300, fps = 10, width = 1000, height = 800)
# anim_save(here::here("faceted_tallest_bin_spectra_wigley.gif"), animated_plot)



# Show the GIF
knitr::include_graphics(here::here("faceted_tallest_bin_spectra_wigley.gif"))

Bodymass spectra abundance peaks

Looking at Individual Years w/ Observable

I like having the ability to look at them one-by-one, but don’t want to make 800 plots or tabs. I’m using this exercise as an excuse to use observable for interactive selections.

# Reformat years as a date - messes with slider
# year_avgs <- mutate(year_avgs, est_year = as.Date(str_c(est_year, "-06-01")))

# Define data to use for js
ojs_define(
  data = wig_l2_aggs %>% arrange(est_year))
  # data = filter(wig_l2_aggs, normalized_abund > 1) %>% 
  #   arrange(est_year))

viewof yearSlider = Inputs.range(
  [1970, 2019], // If not pulling from the data
  {label: "Year:",
   value: 1970, // If not pulling from the data
   step: 1}                    // Step size for numeric slider
)

// Create selections for season and area
//viewof seasonSelect = Inputs.select(["Spring", "Fall"], { label: "Season" })
viewof areaSelect = Inputs.select(["Northeast Shelf", "Gulf of Maine", "Georges Bank", "Southern New England", "Mid-Atlantic Bight"], { label: "Area" })


// // Filtering Function
// filteredData = transpose(data).filter(function(d) {
//   return yearSlider == d.est_year && d.area === areaSelect;
// })
filteredSpring = transpose(data).filter(function(d) {
  return yearSlider == d.est_year && d.area === areaSelect && d.season === "Spring";
})
filteredFall = transpose(data).filter(function(d) {
  return yearSlider == d.est_year && d.area === areaSelect && d.season === "Fall";
})

// Plot with Facets
Plot.plot({
  facet: {
    data: filteredSpring,
    x: d => d.area, // Columns: area
    y: d => d.season  // Rows: season
    //y: { field: "season", domain: ["Spring", "Fall"] }  // Custom row order - doesn't work
  },
  marks: [
    Plot.barY(filteredSpring, { 
    x: "left_lim", 
    y: "normalized_abund",
    fill: "#00608a",
    //stroke: "season",
    stroke: d => d.is_tallest ? "black" : "none", // Change stroke color for highlighted bars
    strokeWidth: d => d.is_tallest ? 4 : 0 // Thicker stroke for highlighted bars
    })
  ],
  x: {
    label: "Body Mass (g)",
    tickFormat: d => `2^${d}` // Format labels as powers of 2
  },
  y: {
   label: "Abundance",
   type: "symlog",
   base: 2
  },
  style: {
    facetPadding: 5
  }
})

// Plot with Facets
Plot.plot({
  facet: {
    data: filteredFall,
    x: d => d.area, // Columns: area
    y: d => d.season  // Rows: season
    //y: { field: "season", domain: ["Spring", "Fall"] }  // Custom row order - doesn't work
  },
  marks: [
    Plot.barY(filteredFall, { 
    x: "left_lim", 
    y: "normalized_abund",
    fill: "#ea4f12",
    //stroke: "season",
    stroke: d => d.is_tallest ? "black" : "none", // Change stroke color for highlighted bars
    strokeWidth: d => d.is_tallest ? 4 : 0 // Thicker stroke for highlighted bars
    })
  ],
  x: {
    label: "Body Mass (g)",
    tickFormat: d => `2^${d}` // Format labels as powers of 2
  },
  y: {
   label: "Abundance",
   type: "symlog",
   base: 2
  },
  style: {
    facetPadding: 5
  }
})

Peak Location Frequencies

Without over-doing the peak locations for multi-modal situations, here are the peak locations and their frequencies for each group:

wig_l2_aggs %>% 
  filter(is_tallest) %>% 
  mutate(est_year = as.numeric(est_year)) %>% 
  ggplot(aes(x = left_lim, fill = season)) +
  geom_histogram(alpha = 0.7, position = "dodge", bins = 9) +
  scale_fill_gmri() +
  scale_x_continuous(
    labels = label_math(expr = 2^.x),
    breaks = c(0:15)) +
    facet_grid(area~season) +
  labs(
    x = "Minimum Body-Size (g)", 
    y = "Year Count",
    title = "Spectra Peak Location Frequency")

Applying the Minimums

From this dataframe we can strip out the key group information and the minimum size as a lookup table. This can be joined and used as a filter as a preparatory step before fitting the group spectra.

# # Make the key
# min_size_key <- wig_l2_aggs  %>%
#   filter(is_tallest) %>%
#   distinct(area, season,  est_year, left_lim, is_tallest) %>%
#   mutate(min_cutoff_g = 2^left_lim) %>%
#   select(-left_lim) %>%
#   mutate(est_year = as.numeric(est_year))
# 
# 
# 
# # Join the data back to the original
# wigley_in_new <- wigley_in %>%
#   mutate(area = "Northeast Shelf") %>%
#   bind_rows(wigley_in) %>%
#   left_join(min_size_key, join_by("est_year", "area", "season")) %>%
#   filter(wmin_g >= min_cutoff_g)
# 
# 
# 
# # Re-Run the MLE method to get new data
# 
# # and again for 16g to 50kg
# peak_chase_results <- group_binspecies_spectra(
#     ss_input       = wigley_in_new,
#     grouping_vars  = c("est_year", "season", "survey_area"),
#     abundance_vals = "numlen_adj",
#     length_vals    = "length_cm",
#     use_weight     = TRUE,
#     isd_xmin       = NULL,
#     global_min     = FALSE,
#     isd_xmax       = NULL,
#     global_max     = FALSE,
#     bin_width      = 1,
#     vdiff          = 2) %>% 
#   mutate(est_year = as.numeric(est_year)) %>% 
#   left_join(area_df)
# 
# 
# shelf_peak_results  <- group_binspecies_spectra(
#     ss_input       = filter(wigley_in_new, area == "Northeast Shelf"),
#     grouping_vars  = c("est_year", "season"),
#     abundance_vals = "numlen_adj",
#     length_vals    = "length_cm",
#     use_weight     = TRUE,
#     isd_xmin       = NULL,
#     global_min     = FALSE,
#     isd_xmax       = NULL,
#     global_max     = FALSE,
#     bin_width      = 1,
#     vdiff          = 2) %>% 
#   mutate(
#     est_year = as.numeric(est_year),
#     area = "Northeast Shelf") %>% 
#   left_join(area_df)
# 
# 
# # Join those together for one file
# moving_peak_spectra <- bind_rows(peak_chase_results, shelf_peak_results)


# # Save them
# write_csv(min_size_key, here::here("Data/processed/wigley_species_l2peaks_key.csv"))
# write_csv(moving_peak_spectra, here::here("Data/processed/wigley_species_l2peaks_bmspectra.csv"))

Plot New Results

# Re-Load them
min_size_key <- read_csv(here::here("Data/processed/wigley_species_l2peaks_key.csv"))
moving_peak_spectra <- read_csv(here::here("Data/processed/wigley_species_l2peaks_bmspectra.csv"))

# Function to get the Predictions
# Flag effect fits that are non-significant  ###
get_preds_and_trendsignif <- function(mod_x){
  modx_preds <- as.data.frame(
    # Model predictions
    ggpredict(
      mod_x, 
      terms = ~ est_year * area * season) ) %>% 
    mutate(
      area = factor(group, levels = area_levels_long_all),
      season = factor(facet, levels = c("Spring", "Fall")))
  
    # Just survey area and year
    modx_emtrend <- emtrends(
      object = mod_x,
      specs =  ~ area*season,
      var = "est_year") %>% 
      as_tibble() %>% 
      mutate(
        zero = 0,
        non_zero = if_else(between(zero, lower.CL, upper.CL), F, T))
  
    # Preds with signif
    modx_preds %>% left_join(select(modx_emtrend, area, season, non_zero))
  
}

# Model linear trends in time
library(nlme)
library(ggeffects)
library(emmeans)
library(tidyquant)

moving_peak_spectra <- moving_peak_spectra %>% 
  mutate(
    area = factor(area, levels = area_levels_long_all),
    est_year = as.numeric(est_year))
spectra_trend_mod <- gls(
  b ~ scale(est_year) * area * season,
  correlation = corAR1(form = ~ est_year | area/season),
  data = moving_peak_spectra)

#check_model(spectra_trend_mod)
# plot(check_collinearity(spectra_trend_mod))


# Get the predictions and flag whether they are significant
bmspectra_trend_predictions <- get_preds_and_trendsignif(spectra_trend_mod) %>% 
  mutate(metric = "bodymass_spectra",
         area = factor(area, levels = area_levels_long_all))



# Make the plot
bmspectra_trend_p <- ggplot() +
   geom_ribbon(
    data = filter(bmspectra_trend_predictions, non_zero == TRUE),
    aes(x, ymin = conf.low, ymax = conf.high, fill = season),
    linewidth = 1, alpha = 0.35) +
  geom_point(
    data = moving_peak_spectra,
    aes(est_year, b, color = season),
    size = 0.8, alpha = 0.5) +
  geom_ma(
    data = moving_peak_spectra,
    aes(est_year, b, color = season, linetype = "5-Year Moving Average"), 
    n = 5, ma_fun = SMA, key_glyph = "timeseries" ) +
  geom_line(
    data = filter(bmspectra_trend_predictions, non_zero == TRUE),
    aes(x, predicted, color = season, linetype = "Significant Trend"),
    linewidth = 1, alpha = 1) +
  scale_fill_gmri() +
  scale_color_gmri() +
  facet_grid(area ~ ., scales = "free") +
  labs(
    title = "Trends in Spectra w/ Shifting wmin",
    x = "Year",
    y = "Exponent of Mass Size Spectra",
    linetype = "Trend",
    color = "Season",
    fill = "Season")
 
# Show plot
bmspectra_trend_p

Things to check in the morning: - autocorrelative structure - are seasons correlated now - what is the relationship of steepness to minimum weght cutoff?

Autocorrelation Check

These are annual chek-ins of the multi-species, multi-cohort community size structure in an area. There should be some memory of the previous six months and/or the previous year.

# Pull out acf stuff
spectra_acf <- moving_peak_spectra %>% 
  unite("regseas", area, season, sep = "XX") %>% 
  split(.$regseas) %>% 
  map_dfr(function(x){
    # drop NA
    x <- arrange(x, est_year) %>% drop_na(b)
    
    # Do ACF
    x_acf <- acf(x$b, plot = F, lag.max = 10)
    
    # Get the signif:
    # https://www.squaregoldfish.co.uk/programming/r_acf_significance.md/
    # Not 100% sure n is the same for ccf as it is for acf, but...
    significance_level <- qnorm((1 + 0.95)/2)/sqrt(sum(!is.na(x$b)))
    
    data.frame(
      "acf"    = x_acf$acf,
      "lag"    = x_acf$lag,
      "sigpos" = significance_level,
      "signeg" = significance_level*-1
    )
    
  }, .id = "regseas") %>% 
  separate("regseas", into = c("area", "season"), sep = "XX")  %>% 
  mutate(
    # Flag when it crosses threshold
    sig_flag = ifelse(acf < signeg | acf > sigpos, T, F),
    # Set Factor Levels
    area = factor(area, levels = area_levels_long_all),
    season = factor(season, levels = c("Spring", "Fall"))) 



# Plot things
ggplot(spectra_acf, aes(lag, acf)) +
  geom_hline(yintercept = 0, color = "gray25", linewidth = 1) +
  geom_vline(xintercept = 0, color = "gray25", linewidth = 1) +
  geom_ribbon(
    aes(ymin = signeg, ymax = sigpos), 
        alpha = 0.2, linetype = 2, fill = "lightgray",
    color = "gray65") +
  geom_col(alpha = 0.65) +
  geom_col(
    data = filter(spectra_acf, sig_flag), color = "black") +
  scale_x_continuous(
    limits = c(-1,10),
    breaks = c(1:9),
    expand = expansion(add = c(0,0))) +
  facet_grid(area~season) +
  labs(
    title = "Spectra Slope ACF",
    x = "Lag (years)", 
    y = "ACF", )

Based on the auto-corrrelation function results we have a similar outlook on the annual autocorrelations.

moving_peak_spectra %>% 
  select(est_year, area, season, b) %>% 
  pivot_wider(names_from = "season", values_from = "b") %>% 
  ggplot(aes(Spring, Fall)) +
  geom_point() +
  geom_smooth(method = "lm") +
  ggpmisc::stat_poly_eq(ggpmisc::use_label(c("P", "R2")), formula = y ~ x) +
  scale_y_continuous(expand = expansion(mult = c(0, 0.5))) +
  facet_wrap(~area, ncol = 2, scales = "free") +
  labs(title = "Spring-Fall Slope Correlations")

Spring and fall communities from the same area remain uncorrelated. This is surprising to me as the noise coming from the presence/absence of cohorts of smaller individuals should be reduced.

moving_peak_spectra %>% 
  ggplot(aes(xmin_actual, b)) +
  geom_point() +
  geom_smooth(method = "lm") +
  ggpmisc::stat_poly_eq(ggpmisc::use_label(c("P", "R2")), formula = y ~ x) +
  scale_y_continuous(expand = expansion(mult = c(0, 0.5))) +
  scale_x_continuous(trans = transform_log(base = 2), labels = label_log(base = 2)) +
  facet_wrap(~area, ncol = 2, scales = "free") +
  labs(title = "Bias from Shifting the Distribution's wmin")

Affinities to Large External Factors

Here I could explore whether the relationships to bottom temperature or landings, but I remain unconvinced this is the right path forward.

# From Andy Beet
landings_annual <- read_csv(here::here("Data/processed/BEET_GARFO_regional_finfish_landings.csv")) %>% 
  rename(area = survey_area,
         est_year = year) %>% 
  left_join(area_df) %>% 
  select(survey_area, area, est_year, total_weight_lb = total_live_lb)

# Get shelf totals
landings_annual <- landings_annual %>% 
  group_by(est_year) %>% 
  summarise(total_weight_lb = sum(total_weight_lb, na.rm = T)) %>%
  mutate(
    survey_area = "Northeast Shelf",
    area = "Northeast Shelf") %>% 
  bind_rows(landings_annual)
  

# Seasonal Bottom Temp
bot_temps <- read_csv(here::here("Data/processed", "trawl_region_seasonal_bottom_temps.csv")) %>% 
  rename(area = survey_area,
         est_year = year) %>% 
  left_join(area_df) %>% 
  select(survey_area, area, est_year, season, bot_temp, overall_avg_btemp, bot_temp_anom)



# Join these
peak_chase_model_df <- left_join(moving_peak_spectra, bot_temps) %>% left_join(landings_annual) %>% 
  mutate(area = factor(area, levels = area_levels_long_all))

# Just drop the random effect part for now
full_model_ar <- nlme::gls(
  # Model
  model = b ~ area * season * scale(bot_temp) + area * scale(log(total_weight_lb)), 
  
  # Data
  data = peak_chase_model_df %>%
    drop_na(total_weight_lb) %>% 
    dplyr::filter(area != "Northeast Shelf"),
  
  # Autocorrelation
  correlation = corAR1(form = ~ est_year | area/season))

# Check the model
library(performance)
# check_model(full_model_ar)


# # Collinearity without interactions
# plot(
#   check_collinearity(
#     nlme::gls(
#         b ~ area + season + scale(bot_temp) +  scale(log_land),
#         data = drop_na(wigley_bmspectra_model_df, total_weight_lb) %>%
#                  filter(area != "Northeast Shelf") %>% 
#           mutate(log_land = log(total_weight_lb)),
#         correlation = corAR1(form = ~ est_year | survey_area/season)
#       )
#     )
#   )


# # confidence intervals on phi
# intervals(full_model_ar)$corStruct

# Clean up the plot:
# Plot the predictions over data
temp_preds <- as.data.frame(
  ggpredict(
    full_model_ar, 
    terms = ~ bot_temp * area * season))   %>% 
  mutate(
    area = factor(group, levels = area_levels_long_all),
    season = factor(facet, levels = c("Spring", "Fall")))



#### Trend Posthoc  ####

# Flag effect fits that are non-significant  ###
temp_emtrend <- emtrends(
  object = full_model_ar,
  specs =  ~ area * season,
  mode = "appx-satterthwaite",
  var = "bot_temp") %>%
  as_tibble() %>%
  mutate(
    zero = 0,
    non_zero = if_else(between(zero, lower.CL, upper.CL), F, T))


# Limit temperature plotting range
actual_range <- peak_chase_model_df %>%
  filter(area != "Northeast Shelf") %>% 
  group_by(season, area) %>%
  summarise(
    min_temp = min(bot_temp)-2,
    max_temp = max(bot_temp)+2,
    .groups = "drop")



# Plot over observed data
local_btemp <- temp_preds %>% 
  left_join(actual_range) %>%
  filter((x < min_temp) == F,
         (x > max_temp) == F) %>%
  left_join(select(temp_emtrend, area, season, non_zero)) %>%
  filter(non_zero) %>%
  ggplot() +
  geom_ribbon(aes(x, ymin = conf.low, ymax = conf.high, group = season), alpha = 0.1) +
  geom_line(
    aes(x, predicted, color = season), 
    linewidth = 1) +
  geom_point(
    data = filter(peak_chase_model_df, area!= "Northeast Shelf"),
    aes(bot_temp, b, color = season),
    alpha = 0.4,
    size = 1) +
  facet_grid(area~., scales = "free", labeller = label_wrap_gen(width = 8)) +
  scale_color_gmri() +
  scale_x_continuous(labels = label_number(suffix = deg_c)) +
  labs(y = "Mass distribution Exponent",
       x = "Local Seasonal Bottom Temperature Anomaly")
local_btemp

# Clean up the plot:
# Plot the predictions over data
f_preds <- as.data.frame(
  ggpredict(
    full_model_ar, 
    terms = ~ total_weight_lb * area * season))   %>% 
  mutate(
    area = factor(group, levels = area_levels_long_all),
    season = factor(facet, levels = c("Spring", "Fall")))



#### Trend Posthoc  ####

# Flag effect fits that are non-significant  ###
f_emtrend <- emtrends(
  object = full_model_ar,
  specs =  ~ area,
  var = "total_weight_lb",
  mode = "appx-satterthwaite"
  ) %>%
  as_tibble() %>%
  mutate(
    zero = 0,
    non_zero = if_else(between(zero, lower.CL, upper.CL), F, T))


# Limit temperature plotting range
actual_range <- peak_chase_model_df %>%
  group_by(season, area) %>%
  summarise(
    min_f = min(total_weight_lb)-2,
    max_f = max(total_weight_lb)+2,
    .groups = "drop")



# Plot over observed data
landings_ar_plot <- f_preds %>% 
  left_join(actual_range) %>%
  filter((x < min_f) == F,
         (x > max_f) == F) %>%
  left_join(select(f_emtrend, area,non_zero)) %>%
  filter(non_zero) %>%
  ggplot() +
  geom_ribbon(aes(x, ymin = conf.low, ymax = conf.high, group = season), alpha = 0.1) +
  geom_line(
    aes(x, predicted, color = season), 
    linewidth = 1) +
  geom_point(
    data = filter(peak_chase_model_df, area != "Northeast Shelf"),
    aes(total_weight_lb, b, color = season),
    alpha = 0.4,
    size = 1) +
  facet_grid(area~., scales = "free", labeller = label_wrap_gen(width = 8)) +
  scale_color_gmri() +
  scale_x_log10(labels = label_log(base = 10), limits = 10^c(0,10)) +
  labs(y = "Exponent of ISD",
       title = "Total Landings (lb.)",
       x = "Local Seasonal Bottom Temperature Anomaly")
landings_ar_plot

Bottom Temperature CCF

Shelf-Wide:

Regionally:

# Function to grab the correlation data and lag data
get_ccf_vector <- function(x,y, lagmax = 10){
  
  # Run the ccf
  ccf_data <- ccf(x, y, plot= F , lag.max = lagmax)
  
  # Get the signif:
  # https://www.squaregoldfish.co.uk/programming/r_acf_significance.md/
  # Not 100% sure n is the same for ccf as it is for acf, but...
  significance_level <- qnorm((1 + 0.95)/2)/sqrt(sum(!is.na(x)))
  
  data.frame(
    "acf"    = ccf_data$acf,
    "lag"    = ccf_data$lag,
    "sigpos" = significance_level,
    "signeg" = significance_level*-1
  )
}


# Run Local CCF on Bottom Temperature
btemp_ccf <- peak_chase_model_df %>% 
  unite(col = "regseas", area, season, sep = "XX") %>% 
  split(.$regseas) %>% 
  map_dfr(~get_ccf_vector(
    x = .x$bot_temp, 
    y = .x$b), .id = "regseas") %>% 
  separate(col = "regseas", into = c("area", "season"), sep = "XX") %>% 
  mutate(
    # Flag when it crosses threshold
    sig_flag = ifelse(acf < signeg | acf > sigpos, T, F),
    # Set Factor Levels
    area = factor(area, levels = area_levels_long_all),
    season = factor(season, levels = c("Spring", "Fall")))


# Plot it
ggplot(btemp_ccf, aes(lag, acf)) +
  geom_hline(yintercept = 0, color = "gray25", linewidth = 1) +
  geom_vline(xintercept = 0, color = "gray25", linewidth = 1) +
  geom_ribbon(
    aes(ymin = signeg, ymax = sigpos), 
        alpha = 0.2, linetype = 2, fill = "lightgray",
    color = "gray65") +
  geom_col(alpha = 0.65) +
  geom_col(
    data = filter(btemp_ccf, sig_flag), color = "black") +
  scale_x_continuous(expand = expansion(add = c(0,0))) +
  facet_grid(area~season) +
  labs(
    title = "Bottom Temperature CCF",
    x = "Bottom Temperature Lag (years)", 
    y = "ACF", )

Landings CCF

Its plausible that landings have a lagged effect so here are the CCFs for that.

# Run Local CCF on Bottom Temperature
landings_ccf <- peak_chase_model_df %>% 
  drop_na(total_weight_lb) %>% 
  unite(col = "regseas", area, season, sep = "XX") %>% 
  split(.$regseas) %>% 
  map_dfr(~get_ccf_vector(
    x = .x$total_weight_lb, 
    y = .x$b), .id = "regseas") %>% 
  separate(col = "regseas", into = c("area", "season"), sep = "XX") %>% 
  mutate(
    # Flag when it crosses threshold
    sig_flag = ifelse(acf < signeg | acf > sigpos, T, F),
    # Set Factor Levels
    area = factor(area, levels = area_levels_long_all),
    season = factor(season, levels = c("Spring", "Fall")))



ggplot(landings_ccf, aes(lag, acf)) +
  geom_hline(yintercept = 0, color = "gray25", linewidth = 1) +
  geom_vline(xintercept = 0, color = "gray25", linewidth = 1) +
  geom_ribbon(
    aes(ymin = signeg, ymax = sigpos), 
        alpha = 0.2, linetype = 2, fill = "lightgray",
    color = "gray65") +
  geom_col(alpha = 0.65) +
  geom_col(
    data = filter(landings_ccf, sig_flag), color = "black") +
  scale_x_continuous(expand = expansion(add = c(0,0))) +
  facet_grid(area~season) +
  labs(
    title = "Total Landings CCF",
    x = "Total Landings Lag (years)", 
    y = "ACF", )